library(tidyverse)
library(tidyNano)
library(themesoo)
library(broom)
library(ggpubr)
library(ggbeeswarm)
library(ggsci)

# devtools::install_github("nguyens7/tidyNano")
# devtools::install_github("soohyuna/themesoo")
my_file <- "PBS_LPS_-ExperimentSummary.csv"


data <- nanoimport(my_file)
## NTA version: 3.2
## Sample name:
## Number of parameters detected: 1
## Dilution factor detected: 1
## Auto name = FALSE
## Custom name: NULL
## Dilution value: 1 (Didn't parse)
data
tidy_df <- data %>%
  nanotidy(sep_var = c("Condition", "Bio_rep","Dilution", "Tech_rep")) %>% 
  mutate(Condition = factor(Condition, levels = c("PBS","LPS")))


tidy_df
tidy_df %>% 
  filter(particle_size < 400) %>% 
  ggplot(aes(x = particle_size, 
             y = True_count, 
             color = Tech_rep)) +
  geom_line() +
  facet_wrap(Bio_rep ~ Condition, ncol = 2)

tech_rep_df <- tidy_df %>% 
  nanolyze(particle_size, Condition, Bio_rep,name = "Tech_rep", param_var = True_count)

tech_rep_df
line_plot <- tech_rep_df %>% 
  ggplot(aes(x = particle_size, y = Tech_rep_mean, color = Bio_rep)) +
  geom_line(size = 1, alpha = 0.7) +
  facet_wrap(~Condition) +
  theme_soo() +
  scale_y_continuous(breaks = seq(0,7E9,1E9)) +
  labs(title = "",
       caption = "(Lines represent biological replicates)",
       x = "Treatment Condition\n",
      y = "Concentration (particles/ml)\n")

line_plot

# ggsave(plot = line_plot, "PBS_LPS_GD16.5_Peripheral_exo_line.png",
#        dpi = 600, width = 7, height = 4, units = "in")
tech_rep_count_df <- tidy_df %>% 
  nanocount(Condition, Bio_rep,Tech_rep,
            name = "Tech_rep", param_var = True_count)

tech_rep_count_df
bio_rep_count_plot <- tech_rep_count_df %>% 
  ggplot(aes(x = Condition, y = Tech_rep_count,fill = Condition, color = Condition)) +
  geom_jitter(width = 0.2, alpha = 0.7, size = 2)  +
  geom_boxplot(fill = NA, width = 0.5, outlier.shape = NA) +
  facet_wrap(~Bio_rep, nrow = 1) +  
  stat_compare_means(data = tech_rep_count_df,
                     method = "t.test",
                     comparisons = list(c("LPS","PBS"))) +
  scale_y_continuous(limits = c(0,7.5E11)) +
  guides(color = FALSE,fill = FALSE) +
  theme_soo() +
  labs(title = "",
       caption = "(Split by biological replicates, points represent technical replicates)",
       x = "Treatment Condition\n",
      y = "Concentration (particles/ml)\n") +
  scale_color_manual(values = c("#0072B5FF","#BC3C29FF")) +
  scale_fill_manual(values = c("#0072B5FF","#BC3C29FF")) 



bio_rep_count_plot

# ggsave(plot = bio_rep_count_plot, "PBS_LPS_GD16.5_Peripheral_exo_bio_reps.png",
#       dpi = 600, width = 8, height = 4, units = "in")
bio_rep_count_df <- tech_rep_count_df %>% 
  nanolyze(Condition,Bio_rep, name = "Bio_rep", param_var = Tech_rep_count)


bio_rep_count_df
count_plot <- bio_rep_count_df %>% 
  ggplot(aes(x = Condition, y = Bio_rep_mean, color = Condition, fill = Condition)) +
  geom_quasirandom(width = 0.1, size = 4) +
  geom_boxplot(alpha = 0.6, width = 0.4, outlier.color = NA) +
  stat_compare_means(data = bio_rep_count_df,
                     method = "t.test",
                     comparisons = list(c("LPS","PBS"))) +
  scale_y_continuous(limits = c(0,6.5E11)) +
  guides(color = FALSE, fill = FALSE) +
  theme_soo() +
  labs(title = "",
       x = "Treatment Condition",
      y = "Concentration (particles/ml)\n") +
  scale_color_manual(values = c("#0072B5FF","#BC3C29FF")) +
  scale_fill_manual(values = c("#0072B5FF","#BC3C29FF")) 


count_plot

mypal <- pal_nejm("default")(2)
mypal
## [1] "#BC3C29FF" "#0072B5FF"
# ggsave(plot = count_plot, "PBS_LPS_GD16.5_Peripheral_exo.png",
#        dpi = 600, width = 6, height = 5, units = "in")
count_plot %>% 
  plotly::ggplotly()
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomSignif() has yet to be implemented in plotly.
##   If you'd like to see this geom implemented,
##   Please open an issue with your example code at
##   https://github.com/ropensci/plotly/issues
shapiro.test(bio_rep_count_df$Bio_rep_mean) %>% 
  tidy()
t.test(Bio_rep_mean~Condition,data = bio_rep_count_df) %>% 
  tidy()